BiocManager::install('SNPRelate', Ncpus = 6)
library(SNPRelate)

# basic index & read -------
snpgdsExampleFileName() |>
  snpgdsSummary()

genofile <- snpgdsExampleFileName()|>
  snpgdsOpen()

genofile |>
  index.gdsn('sample.id') |>
  read.gdsn() |>
  head()

# Take out genotype data for the first 3 samples and the first 5 SNPs
genofile |>
  index.gdsn('genotype') |>
  read.gdsn(start = c(1,1), count = c(5,3))

genofile |>
  index.gdsn('snp.rs.id') |>
  read.gdsn() |>
  head()

# Read population information
genofile |>
  index.gdsn(path="sample.annot/pop.group") |>
  read.gdsn() |>
  table()

# Close the GDS file
snpgdsClose(genofile)

# convert PLINK file to GDS -------
# The PLINK BED file, using the example in the SNPRelate package
bed.fn <- system.file("extdata", "plinkhapmap.bed.gz", package="SNPRelate")
fam.fn <- system.file("extdata", "plinkhapmap.fam.gz", package="SNPRelate")
bim.fn <- system.file("extdata", "plinkhapmap.bim.gz", package="SNPRelate")

snpgdsBED2GDS(bed.fn, fam.fn, bim.fn, "test.gds")

snpgdsSummary("test.gds")

# convert VCF to GDS ---------
## native method ------
# The VCF file, using the example in the SNPRelate package
vcf.fn <- system.file("extdata", "sequence.vcf", package="SNPRelate")

# Reformat
snpgdsVCF2GDS(vcf.fn, "test.gds", method="biallelic.only")

snpgdsSummary("test.gds")

## SeqArray method ---------
# It is strongly suggested to use `SeqArray` for large-scale WGS variant data.
library(SeqArray)

# the VCF file, using the example in the SeqArray package
vcf.fn <- seqExampleFileName("vcf")

# convert, save in "tmp.gds" with the default lzma compression algorithm
seqVCF2GDS(vcf.fn, "test.gds")

# It is suggested to use seqGetData() to take out data from the SeqArray file
# since this function can take care of variable-length data and multi-allelic genotypes, 
# although users could also use read.gdsn()
genofile <- seqOpen("test.gds")

genofile

# take out sample id
samp.id <- genofile |>
  seqGetData('sample.id') |>
  head()

samp.id

# take out variant id
variant.id <- genofile |>
  seqGetData('variant.id') |>
  head()

variant.id

# get "chromosome"
genofile |>
  seqGetData('chromosome') |>
  table()

# get "allele"
genofile |>
  seqGetData('allele') |>
  table()

# get "annotation/info/GP"
genofile |>
  seqGetData('annotation/info/GP') |>
  head()

# set sample and variant filters
genofile |>
  seqSetFilter(sample.id=samp.id[c(2,4,6)])

# or
genofile |>
  seqSetFilter(sample.sel = c(2,4,6))

genofile |>
  seqSetFilter(variant.sel = 1:4)

# get "allele"
seqGetData(genofile, "allele")
## [1] "T,A" "G,A" "G,C" "A,G"

# get genotypic data
# 3-dimensional data: variant, sample, allele
seqGetData(genofile, "genotype")

# get the dosage of reference allele
seqGetData(genofile, "$dosage")
##       variant
## sample [,1] [,2] [,3] [,4]
##   [1,]    2    1    2    2
##   [2,]    2    2    2    2
##   [3,]    2    2    2    2

# close the file
seqClose(genofile)

# Calculate Fst value ------------------
# Get sample id
sample.id <- genofile |>
  index.gdsn("sample.id") |>
  read.gdsn()

# Get population information
#   or pop_code <- scan("pop.txt", what=character())
#   if it is stored in a text file "pop.txt"
pop_code <- genofile |>
  index.gdsn(path="sample.annot/pop.group") |>
  read.gdsn() 

# Two populations: HCB and JPT
flag <- pop_code %in% c("HCB", "JPT")
samp.sel <- sample.id[flag]
pop.sel <- pop_code[flag]
v <- genofile |>
  snpgdsFst(sample.id=samp.sel, population=as.factor(pop.sel),
               method="W&C84")

v$Fst        # Weir and Cockerham weighted Fst estimate
v$MeanFst    # Weir and Cockerham mean Fst estimate
summary(v$FstSNP)

# Multiple populations: CEU HCB JPT YRI
#   we should remove offsprings
father <- read.gdsn(index.gdsn(genofile, "sample.annot/father.id"))
mother <- read.gdsn(index.gdsn(genofile, "sample.annot/mother.id"))
flag <- (father=="") & (mother=="")
samp.sel <- sample.id[flag]
pop.sel <- pop_code[flag]
v <- snpgdsFst(genofile, sample.id=samp.sel, population=as.factor(pop.sel),
               method="W&C84")